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I . INTRODUCTION 


This report constitutes a final report on research into the effects 
of hypersonic nonlinear aerodynamic loading on panel flutter, supported 
by the National Aeronautics and Space Administration under Grant NGR 05 - 
020-102, from September 1, 19^5 to August ^1, I 9 TO . The technical moni- 
tor of this grant was Mr. Peter A. Gaspers, of the Nonsteady Phenomena 
Branch at Ames Research Center, whose support and counsel are gratefully 
acknowledged , 

The initial motivation for this investigation was provided by the 
results of some high -Mach -number panel-flutter experiments described in 
Ref, 1, In these experiments, in-plane tension was used to stabilize a 
panel until steady-state tunnel conditions were reached. The tension was 
reduced until flutter occurred and then increased until flutter ceased. 

It was noticed that the tension at which flutter ceased was in many cases 
consistently higher than that for which flutter occurred initially, so 
that there was an effect analogous to hysteresis evident in the relation- 
ship between flutter speed and in-plane tension. Such an effect is con- 
sistent with the effects of nonlinear aerodynamic loading at hypersonic 
speeds, which are easily and accurately represented for panels by non- 
linear piston-theory aerodynamics. In general, these nonlinear aero- 
dynamic loads are of the ^‘soft" type, in that they produce forces aug- 
menting any transverse panel motion. On the other hand, the more fa- 
miliar geometric panel, nonlinear effects are of the "hard" type, in that 
they produce in-plane panel loads that resist transverse displacement. 

It thus seemed conceivable that there would be circumstances, perhaps 
exemplified by the experiments in Ref, 1, where the aerodynamic non- 
linear effects would in part be dominant and lead to the hysteresis 
effect observed. There is also the broader question of stability to 
large disturbances; Refs. 2 and 3 deal with this in different ways. 

In Ref. 2, Bolotin and co-authors showed that nonlinear hypersonic 
aerodynamic loads can produce situations where a panel is unstable in 
a parameter region where stability would be predicted with linear aero- 
dynamic loading, provided the initial disturbance is of sufficient 
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magnitude. In. Ref. 3^ Librescu was able to demonstrate a similar effect, 
albeit not so precisely, with an analysis based on Lyapunov stability 
theory. This question was pursued further, with more realistic panel 
parameters, as part of the study described in this report. The effects 
of nonlinear aerodynamic loads on panel postcritical response were also 
surveyed, and results have been published in Ref. h. The effects of 
aerodynamic nonlinearities on stability were evaluated by determining 
constant-initial-energy amplitude -sensitive stability boundaries and 
comparing them with the corresponding linear stability boundaries. Pre- 
liminary results were presented in Ref, H, and a parametric survey was 
presented in Ref. 5 . A final section treats an attempt to develop an 
alternative method of analysis for systems where amplitude -sensitive 
instability is possible; much of this section is the work of Dr, Sayed 
D. Hassan. More detailed summaries of work performed during the first 
four years of the grant period will be found in Refs. 7-10. 

Another important effort supported by the grant concerned the ef- 
fects of a viscous boundary layer on unsteady panel aerodynamic loading. 
This work is described in detail in the Ph.D. dissertation of Dr. J. I. 
Lerner, and is incorporated into this report by reference (Ref. 6). 

The author is greatly indebted to Professors Holt Ashley, jian Mayers, 
and Krishnamurty Karamcheti for many helpful discussions, to Mr. James 
Stein for valuable computer-programming assistance, and to Dr. Ahmed A. 
Kamel for consultations concerning his perturbation method. 
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II. EQUATIONS OF MOTION 


Consider the two-dimensional panels or plate -column, illustrated in 
Fig, 1. The supports are hinged, and the in-plane motion at one end is 
resisted by a distributed spring of constant K .* The panel is loaded 
by a static pressure difference Ap , an unsteady pressure difference 
p(x, t) - p^ , an initial static in-plane applied load , and a sub- 
sequent time -varying in-plane applied load • The unsteady pres- 

sure is approximated in hypersonic flow by a third-order piston-theory 

- - 

M . M ) 


expression 


P - P. 


M 


1 ^ 
U St 




(7 + 1 )M 

V 


+ 


(7 + 1)M^ 

12 



( 2 . 1 ) 


A third-order term is included, because such a term will result in work 
done through a simple-harmonic cycle of panel motion, whereas the second- 
order terms will not. The pressure is assumed to act normal to the in- 
stantaneous panel surface. Since the structural representation of the 
panel allows for significant rotation of panel elements about a span- 
wise axis, order-of-magnitude consistency then dictates including as 
well as transverse aerodynamic - loading an in-plane aerodynamic load 
[(P“Poo)l “ , where the subscript L denotes the linear 

portion in Eq. (2,1). 

The panel transverse displacement is represented as a series of 
assumed modes satisfying the geometric boundary conditions of zero 


See the Nomenclature for definitions of symbols. 
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displacement and curvature at each end: 


N 

w(x,t) = ^ 

k=l 

A consistent assumed-mode expression 
given by 

u(x;,t) := [h^ + b^Ct)] 

Here b is the initial panel in-plane displacement at x = a due to 
R 

the application of the in-plane load R j and b (t) is the in-plane 

X vy 

displacement at x ^ a resulting from the subsequent unsteady panel mo- 
tion. The reasons for this division will be discussed below, 

Hamilton's principle is now used to derive the Euler-Lagrange equa- 
tions of motion, as outlined in Ref. 7. The potential energy of the sys- 
tem is given by 

^w 
^x 


( 2 .!*) 

This expression represents the strain energy of the panel diminished by 
the potential of any conservative external loads. The structural repre- 
sentation of the panel is consistent with the assumptions of small strain, 
a linear stress -strain law, and element rotations whose squares are small 
relative to unity (see Ref. U), and it leads to the familiar von Karman 
structural operator. In the potential, it has been assumed that the re- 
straint spring is not attached until after the load R is applied; this 

X 

is an artifice for ensuring that in-plane tension is in the panel even 






a^(t) sin ^ (2.2) 

for the in-plane displacement is 


2N 

k=l 


h 



when K = “ , The unknown b is needed simply to account for the 

corresponding panel displacement. This formulation also implies that 
the panel is initially restrained from buckling, should supercritical 
compressive in-plane loading be applied. 

The kinetic energy is given simply by 

0 

The kinetic energy associated with both in-plane motion and rotary in- 
ertia has been neglected. In-plane inertia would have to be included 
if parametric stability of the panel, under the load ^Q(b) j were to 
be studied. Here is included only to model the varying tension 

used in the experiments of Ref. 1, where in-plane inertia clearly was 
not a factor. Rotary-inertia effects are important only if the wave- 
length of the panel flexural mode is comparable in magnitude to the 
panel thickness; such small wavelengths are not expected for panel 
flutter. 

The transverse generalized forces are 

a 

" f t- (p “Poo) + Ap] sin ^ dx , k= l,2, ...,N (2.6) 

0 

Note that unsteady cavity effects, which in some cases are quite signif- 
icant, are not taken into account here. The in-plane generalized forces 
are 

a 

= / t(p-P«.)L- 

0 

a 

= / '(P-P»>L - 

0 


Ap ] ^ — dx 
a 


(2.7) 


^ . kTTX 

sxn 

^x a 


sin dx , k 


= 1,2/...,2N 


(2.5) 
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The assumed-mode series for u and w are inserted in the expres- 
sions for potential energy^ kinetic energy, and generalized forces. The 
Euler -Lagrange differential equations then become, after considerable 
manipulation. 
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Multiple summations in these equations are indicated by multiple indices 
on a single summation sign^ and the terms are arranged so that any com- 
bination of indices giving a zero denominator is to be excluded from the 
summation, (This convention will be observed throughout this report.) 
Equations (2.B) and (2,9) ^^re used to eliminate 

in Eqs . (2,11)_, and this produces a set of second-order, quasilinear 
ordinary differential equations governing the a^ . After nondimension- 
alization and copious manipulation, this set becomes: 
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With the exception of a slight difference in the definition of a , and 
the time-varying in-plane load , the terms up to and including the 

first appearance of are the same as those derived by Dowell (Ref. 11): 

the von Kirmdn nonlinear structural operator for the panel_, with linear 
pis ton-theory aerodynamic loading. The system parameters for this com- 
bination of terms are \ , a , ^ ^ and R^ Furthermore, 

Dowell pointed out that these equations can be recast so as to eliminate 
the explicit appearance of o: , so that only one (nonzero) value of Q! 

needs to be considered. With the addition of the nonlinear aerodynamic 
terms, the situation becomes more complicated, A new parameter, h/a , 
appears explicitly, and p. and M must be specified separately. Also, 
the explicit dependence on a can no longer be eliminated. The nonlinear 
aerodynamic terms appear in Eqs , (2,12) in the following order: first, 

those arising from the terms in (F^)^ dependent on (^w/^x)^ , (^w/^x) , 

, and (^w/^t) , respectively; secondly, those arising from the 

terms^n ^F^)^ and (F^)^ dependent on (^w/^x)^ , (^w/^x) , (^w/^t) , 
and Ap (^^) , respectively; finally, those arising from the terms in 

dependent on ($w/^x)^ , Note also that the panel geometric non- 

linear terms are cubic in the a^ , whereas the nonlinear terms from 
linear in h/a and quadratic in the a^ , and those from 
^^x^O ^^x^k quadratic in h/a and cubic in the 

For the- complete equations, the functional dependence of the panel 
transverse displacement can be stated as 

w(x,t) = hm T ; Aq , 1 , n, M, a, 7 , Ap , (2,X5) 

The method of solution is to integrate Eqs. (2.12) from given initial con- 
ditions and to observe the resultant panel motion versus time. The 

use of piston theory, with its point-function relation between pressure 
^nd panel motion, makes this possible. The initial state-space vector 
Aq will not. always govern the asymptotic behavior of w ; instances 
where this does occur will be discussed below. The computer program to 
integrate these equations was set up so as to allow individual nonlinear 
aerodynamic terms to be left out or included, as desired. 
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III. VERIFICATION; PRELIMINARY ASSESSMENT 
OF AERODYNAMIC NONLINEARITIES 

I ' ■ 

5.1. Comparison with Previous Results 

j 

In order to check the numerical integration, a few test cases were 
run with linear aerodynamic loading, and the results were compared with 
those from Ref , 9» Figures 2 and 5 show the dimensionless panel dis- 
placement at the three-quarter chord vs. dimensionless time for values 
of system parameters noted on the figures. Figure 2 corresponds di- 
rectly to Fig, 2 of Ref. 11, and the limit-cycle amplitude taken from 
Fig. 5 checks that given in Fig. 5 of Ref. 11. All the other test cases 
showed similar good agreement, 

5,2 Behavior with Zero System Damping ^ ; 

Some interesting and unanticipated problems were encountered when 
the equations with linear aerodynamic loading were solved for zero sys- 
tem damping. Figures and 5 show the panel response for the same set 
of initial conditions and system parameters used for Figs, 2 and 5, 
except that p = 0 . After an initial transient shown in Fig, the 

response becomes periodic but not simple harmonic, as is evidenced in 
Fig. 5 , The response curve illustrated in Fig. 5 continues indefinitely 
without change, and the peak amplitude is 1 different from that shown in 
Fig. 5. That the system does demonstrate instability for this set of 
parameters is illustrated in Fig. 6, Here the initial displacement at 

I * 

the three-quarter chord is smaller than that in Figs. 2-5 by a factor 
of approximately four, and the amplitude clearly grows with time. It 
then decreases, and the same process is repeated again (but is not shown 
here). Note also that the maximum amplitudes are different for the dif- 
ferent initial conditions . This does not mean, however, that a unique 
limit-cycle amplitude and frequency are not associated with a given set 
of (supercritical) system parameters. In Fig. the panel response is 
shown for N = 2 and for supercritical system parameters ., The initial 
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conditions in this case were obtained from a harmonic-balance solution. 

It is seen that there is no initial transient, and the panel response 
continues at the same amplitude. Solutions for other values of system 
parameters exhibit the same behavior, and it can therefore be concluded 
that the method of calculating the panel motion with time will produce 
a limit cycle of constant amplitude for zero system damping only if the 
initial conditions correspond exactly to the limit -cycle modal ampli- 
tudes, For any other initial conditions, the panel oscillates between 
stable and unstable states. If, on the other hand, the system parameters 
are subcritical, the panel will oscillate without decaying, and the peak 
amplitude is determined by the initial conditions. Figures 8 and 9 show 
the panel displacement at the three-quarter chord for a subcritical value 
of k and different initial conditions. Here it will be observed that 
the peak amplitudes never exceed the initial amplitudes, although the mo- 
tion is not simple harmonic. 

These results can also be interpreted with the aid of some general 
stability considerations for autonomous systems (Ref. 12). The panel 
equations for zero system damping describe a system whose state is 
uniquely determined by the 2N modal ainplitudes and velocities a, , 

iC 

a, , k = 1,2, ...,N , or in other words by a point in the 2N-dimensional 

^ 2N . . 

space E . The origin of this space is clearly an equilibrium point 

3 

— it corresponds to the panel in a flat, undisturbed state — and it is 
desired to examine the stability of the panel in the neighborhood of this 

state. In brief, the origin is stable if for any sphere S (r) of radius 

2N ' 

R in E , centered on the origin, there exists another sphere S(r) 
of radius r g R such that any motion originating in S(r) remains in 
S(r) ever after. This stable behavior is exhibited in Figs. 8 and 9 , 
Reducing the initial amplitude reduces the maximum resultant amplitude, 
so for any given bound on the amplitude the initial amplitude can be re- 
duced to keep the resultant motion within the bound. On the other hand, 
if such a sphere S(r) cannot be found, the origin is unstable. Figures 
5 , and 6 . illustrate this unstable situation. It is evident that there 
is an amplitude that the resultant motion will exceed, no matter how small 
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(but finite) the initial amplitude is. These same conclusions can be 
drawn when the initial conditions are broadened to include nonzero in- 
itial velocities, and the reader is referred to Ref. 12 for the full, 
precise definitions of stability and instability. 

3.5 Isolated Effect of Aerodynamic Nonlinearities 

The effects of aerodynamic nonlinearities were first considered by 

introducing into the equations of motion only the term from (F ), pro- 

— 2 ^ ^ 
portional to (^w/^x) , with the in-plane restraint parameter o: set 

to zero, so that no panel geometric nonlinear terms were present. For 
this single nonlinear aerodynamic term, only the new system parameter 
M(h/a) is required. In Fig, 10, the panel is observed to be stable 
for X = 330 and initial conditions given by “ 1*71 . 

In Fig. 11, the panel diverges with all system parameters unchanged and 
a slight increase in the initial deformation, to a^(0) = - 22(0) = 1,72 , 
The divergence in this case is into the cavity, as expected, since an 
aerodynamic pressure proportional to (^w/^x) produces an increased 
pressure on the free-streara side of the panel for any panel displacement. 
Here, then, the asymptotic panel behavior is clearly governed by the in- 
itial conditions. On the other hand, with linear aerodynamic theory the 
panel would eventually return to its initial flat state, no matter what 
initial conditions would be given. (The classical linear critical value 
of X for this case is approximately 5^5 “ see Ref. 11, for example.) 
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IV. EFFECT OF AERODYNAMIC NONLINEARITIES 


ON POSTCRITICAL RESPONSE 


Of interest here is how the aerodynamic nonlinear terms affect the 
motion and peak stress levels of the panel as it oscillates at a super- 
critical value of X , Comparisons were made between the limit-cycle 
frequency and amplitude obtained with the nonlinear aerodynamic terms 
and those obtained with linear aerodynamic terms. Various combinations 
of system parameters roughly characteristic of those found in practice, 
were tried. Various combinations of nonlinear aerodynamic terms were 
also tried, in an effort to identify those terms that would be important. 
It was soon found that for the system parameters surveyed only two had a 
significant effect — those resulting from terms in (F ) proportional 
to (^w/^x) and (^w/^x) (^w/^t) , With only these aerodynamic non- 
linear terms in Eqs . (2.12), further simplification can be achieved 
through the transformation used by Dowell in Ref. 11. With = (a)^ , 

^0 "" ^0 ^ ^ ~ ^ parameter a can be eliminated 

explicitly from the equations, and in addition to Dowell's parameters a 
single new one appears — Mh/a(o:)^ . This parameter then serves to mea- 

sure the relative importance of nonlinear aerodynamic effects, since it 
ratios the principal factor governing the nonlinear aerodynamic terms to 
the one governing the panel geometric nonlinear terras. The functional 
dependence of the panel transverse displacement simplifies to 


■ . 1 , ■ 1 . 

w(x,t) =, [h/(a)2] Fn[x/a , T ; ^ , X , h/m, r , Mh/a(a)^, ^ , R^, R ] 

(•*. 1 ) 

Postcritical panel response with linear aerodynamic loading was com- 
pared to that with nonlinear aerodynamic loading for various combinations 
of system parameters, A maximum value of 0.05 was taken for Mh/a (say, 
h/a = 0.005 and M == 10.0 ), A typical comparison is illustrated in 
which compares the growth in peak limit— cycle amplitudes at 
x/a = O.T 5 as a function of X for linear and nonlinear aerodynamic 
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loading. Values of other system parameters are given in the caption. 

For the linear aerodynamic case, the panel oscillates as far into the 
airstream as it does into the cavity. With nonlinear aerodynamic loads, 
the peaks into the cavity (negative w) are greater, and those into the 
flow (positive w) are less. This effect is caused primarily by the non- 
linear aerodynamic term dependent upon (^w/^x) , which provides an 

overpressure, tending to push the panel into the cavity, as a result of 
any deviation of the panel from its flat initial position. However, any 
changes brought about by the nonlinear aerodynamic loading are quite 
small, even well into the supercritical regime, and the frequency of 
oscillation is virtually unchanged. 

Since the in-plane restraint parameter a is the key parameter 
governing the influence of the panel geometric nonlinearity, one might 
expect that greater differences would be observed as OL is reduced. A 
linear -nonlinear aerodynamic comparison is shown in Fig. I5 for X ^ 55O 
and varying a . It can readily be observed that the nonlinear aero- 
dynamic terras do have a relatively greater effect, to the extent that the 
peak displacement into the external flow at x/a = 0.75 fs reduced some 
25^ for a =0,1 . On the other hand, the corresponding peak displace- 

ment into the panel is increased by very little, and the frequency re- 
mains virtually unaffected. 

Another important consideration is the change in stress caused by 
nonlinear aerodynamic loading. The stress in the panel can be written 
in terms of modal amplitudes and other system parameters as (Ref. 11) 


.a 


1 -V 
E(h/a) 


N 


2 x 


- V 


k7IX 


a, srn 
k 


k=l 


R + (1 -a)R^ a 
+ ^ 


N 


12 




(4.2) 


k=l 


The maximum or minimum of cr at any instant of time occurs for 

A 

z = ± h/2 , so it is seen that the stress distributions for maximum 
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and minimum stress will plot as curves symmetric about a mean, the stress 

due to stretching, given by the constant terms in Eq. (*<.2). Figure l^i 

compares these stress distributions in the panel, at the instant the 
displacement at x/a - 0,75 reaches a peak, for a = 0.1 and X = 550 . 
The maximum tensile and compressive stresses occur at x/a ^ 0,85 i 
the nonlinear aerodynamic loading increases these stresses by at most 5^« 

For aluminum, E 10 psi and v = I/3 > with h/a = 0,005 y ^ stress 

ff of 60 corresponds to a dimensional stress cr of. 16,900 psi, which 

X X 

is well below the yield stresses of approximately 60,000 psi in tension 
and ^0,000 psi in compression. Convergence studies on stress by Dowell 
(Ref. 9) indicate that six modes (N = 6) may not be adequate for 
stresses, although this number is satisfactory for displacements. How- 
ever, it is not felt that using enough modes to converge the stresses 
would alter the conclusion that the assumption of stress levels within 
the linear range is a valid one. 

Another key assumption in the derivation of the panel equations of 
motion is that (^’w/^x) is everywhere much less than unity. In terms 
of modal amplitudes, this is 



This distribution was calculated for various large -amplitude panel states 
Generally, the largest values were at x = a , and the largest of these 
calculated was approximately O.Ol, which was calculated from the mode 
shape corresponding to the nonlinear-aerodynamic stress distribution of 
Fig, 1^, Values calculated for other cases and over other portions of 
the panel were much smaller, thereby strongly suggesting that the 
’'moderate-rotation'* assumption is also justified. 

One final check involved testing the accuracy o.f the numerical in- 
tegration. This was done by integrating the equations of motion backward 
in time, with initial conditions given by the state of the panel at some 
instant during a previous calculation. Panel motions with time could 
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then be compared. The parameter in the numerical -integration subroutine 
governing acceptable relative error between integration steps was kept 
small enough so that no appreciable differences arose between forward 
and backward integrations over the longest time intervals contemplated. 

The results presented in this section are representative of the re- 
sponses calculated for a variety of conditions. From this evidence, it 
can be concluded that the influence on postcritical response of nonlinear 
aerodynamic loading, of the type considered herein, is minor. The post- 
critical motion is generally such that the aerodynamic terms quadratic 
in ^w/^x do very little net work over a cycle, and the cubic terms, 
which would do work over a cycle, are not important unless Mh/a is 
unrealistically large. It was observed, however, that transient mo- 
tions — such as those that were calculated before steady-state ampli- 
tude and frequency were attained — were affected markedly by the non- 
linear aerodynamic terms. It is therefore entirely possible that these 
terms would be very important when panel transient response is studied 
at hypersonic speeds. The two nonlinear aerodynamic terms used in this 
section were the only ones retained for subsequent studies, except for 
occasional check runs to ensure that other terms were not contributing 
significantly to the results. 
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V. EFFECT OF AERODYNAMIC NONLINEARITIES 
ON STABILITY 


5,1 Comparison with Experiment 

As a first step towards reproducing the experimental conditions of 
Ref. If as discussed in Sec. I, it was. decided to study the panel re- 
sponse as the in-plane applied load was varied across the linear sta- 
bility boundary. This involves scheduling the time-varying load 
with time; a typical case is shown in Fig. I5 . The panel was initially 
set in motion at point A , on the stable side of the linear stability 
boundary. Then the in-plane load was decreased to point B , on the 
unstable side, and held there until the motion of the panel was estab- 
lished, Finally, the load was increased to point A , where the non- 
linear aerodynamic terms dictate the resultant stability characteristics. 
Figure I6 presents a time history of the displacement at x/a = 0.75 for 
values of X and in-plane load corresponding to point A in Fig. I5. 

The motion shown is that which results after the in-plane load has been 
cycled from point A to point B and back to point A . Figure I6 
shows that the energy imparted to the panel while the load is at point 
B is enough to cause instability with the load at point A . To 
stabilize the panel, it would be necessary to increase the load even 
more . 

The scenario described above is consistent with the experimental 
observations of Ref. 1, but the parameters used were not comparable. 

The next step was then to attempt a closer comparison by looking for an 
amplitude -sensitive instability with parameters corresponding to the ex- 
perimental conditions. The following parameters were chosen: R = I60 , 

X = 2000 , M = 10 , }i = 0.1 , = 0 , and h/a = 0.0005^ . These values 

give a condition that is just on the stable side of the linear stability 
boundary for a panel on hinged supports. Experimental edge conditions 
would be better represented with clamped supports, but the theoretical 
differences between stability boundaries for these two edge conditions 


22 



are not really significant when the in-plane applied loads are as large 
as they were in the experiments (see Ref. l). The remaining unknown 
parameter is OL . Various initial amplitudes were used, appropriate 
to flutter amplitudes observed experimentally (w 10 in some cases), 
and a was varied in order to see if these initial amplitudes would 
produce an unstable panel motion. The only unstable motion that could 
be produced was an oscillatory but divergent one, for values of a on 
the order of 10 Hiese results indicate that nonlinear aerodynamic 

influences are not the cause of the experimentally observed behavior, 
since the value of o: needed and the corresponding calculated unstable 

motion are not consistent with the experimental setup or observations . 

5,2 Constant-Initial-Energy Stability Boundary 

After it was demonstrated that energy levels capable of causing 
amplitude-sensitive instability could be generated by unstable panel 
motion near the linear stability boundary, it was decided to consider 
how the linear stability boundary presented, say, in the ^ ~ 
plane with other parameters fixed — would be changed for a given level 
of initial excitation. 

Portraying analytically the dependence of the amplitude-sensitive 
— > 

instability on would be a formidable, if not impossible, task, 

since it has 2 N elements. One simplified approach is suggested by the 

work of Dimantha and Roorda (Ref, 15 )« They investigate the stability 

of nonlinear nonconservative systems with the direct method of Liapunov, 

with Zubov's procedure for constructing the Liapunov functional. It is 

not at all clear whether this method is applicable to the panel stability 

problem; however, the ideas in Ref, I3 do at least suggest a meaningful 

procedure for the numerical experimentation described in this report. 

In principle, one can determine a stability boundary for the panel in 

— ^ 

terms of the elements of A^ for fixed system parameters. This boundary 
can be viewed as a hypersurface S in the phase space determined by the 
elements of A^ , (Note that the origin of this phase space represents 
the panel in its flat undisturbed equilibrium position,) Any combination 
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of initial conditions that plots on one side of S will result in un- 
stable panel motion, while a combination that plots on the other side 
will result in stable motion, Dimantha and Roorda proposed to calculate 
the minimum total panel energy on S such that this minimum value 
will determine a hypersphere that just touches S at one or more points 
and is everywhere else on the stable side of S . The energy then 

provides an upper bound for the initial disturbance energy such that the 
resultant panel motion is stable. In similar fashion, Dimantha and 
Roorda also proposed to determine the maximum total energy on S , 
thereby determining a hypersphere of energy level that gives a 

lower bound on the initial disturbance energy for unstable motion. The 
reader is referred to Ref. I 5 for the full details of this analysis. It 
suffices to say here that these ideas suggest deteirmihing a nonlinear 
stability boundary in the X *- R^ plane for constant initial energy 
and comparing this boundary with the linear stability boundary, which is 
independent of the initial conditions. 

The dimensionless panel energy is given by (Ref. lO): 
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(5.1) 


Note that the time-varying applied load, R^ , is here set to zero. 

The initial choice of energy level was that corresponding to supercritical 
panel motion near the linear X - R stability boundary. It was further 
decided, purely arbitrarily, to set a (o) = - a (O) and to let these be 
the only nonzero components of . With = 1?50 and a = 0.1 , 

the variation of (o) = - a (o) with R is given in Fig, 17. 

The remaining system parameters were then chosen, and a new stability 
boundary was determined, as is shown in Fig, I8. This boundary was ob- 
tained by integrating the panel equations of motion, with initial 
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conditions determined as discussed above^ and observing whether or not 
the calculated panel motion persisted or died out past the initial 
transient. The complicated nature of the dependence of this stability 
boundary on the initial energy (and, ultimately, on ) can be in- 

ferred from the shape of the boundary. As noted in Ref, 5^ the onset 
of the instability over virtually the whole boundary is characterized 
by a strong traveling -wave component, quite similar to the ’’periodic 
but non- simple -harmonic" behavior noted by Dowell (Ref, 11) for large 
compressive values of R^ with linear aerodynamic loading. Observe 
also that the maximum reduction in stability occurs for zero or slightly 
positive (tensile) values of R . The behavior of the nonlinear sta- 
bility boundary near R = - 0.5 77 has not been explained. 

5 .5 Parametric Survey 

The procedure of the previous subsection was then used to determine 

unstable regions with variations in other system parameters « Figure I 9 

gives the variation of the value of X where instability first occurs 

with both different energy levels and different modal content. For 

E = IT 50 ' y there are only minor differences in the critical value of 

— > 

X among the three sets of nonzero elements of A^ that were examined. 
Almost tripling the initial energy for one of these sets produced little 
change in the critical value of X , For = 55^0 , there is an un- 

stable band bounded by stable regions, suggesting the same sort of con- 
tortions in the stability boundary near R^ = 0 as was observed in 
Fig. 18 near R = - 0.5 77^ . 

Figure 20 presents the effects of both positive and negative static 
pressure differences for R^ = 0 and the initial conditions of Fig, I 8 . 
Negative values of produce a static pressure tending to push the 

panel into the cavity, reinforcing the deformation caused by the aero- 
dynamic term proportional to (^w/^x) and slightly adding to its de- 
stabilizing influence. Just the opposite occurs for positive values 
of . 
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As would be expected^ the parameter Mh/a(a)^ has a strong in- 
fluence on the amplitude-sensitive instability, as is shown in Fig. 21. 
Increasing values of Mh/a(a)^ , which represent increased relative 
importance of the nonlinear aerodynamic terras, result in monotonically 
decreasing critical values of X . The slopes of these curves at the 
X axis give a measure of the importance of nonlinear aerodynamic effects. 
Generally, they should all be zero for zero initial energy and decrease 
as the initial energy increases, but the results in Fig. I 9 suggest a 
complicated picture, and any unsupported generalization must be viewed 
with caution. 



VI. A NEW METHOD OF ANALYSIS 


6,1 Kamel ’s Perturbation Method 


Although it produces a great amount of information^ direct integra- 
tion of the equations of motion has certain disadvantages. The informa- 
tion produced is essentially similar to that produced experimentally, and 
the analyst is forced to conduct a number of "numerical experiments" in 
order to determine the influence of the system parameters on the panel's 
behavior. Clearly, it would be extremely worthwhile to have a method of 
analysis that would permit the direct analytical determination of the 
information desired. For example, such a method would be very useful in 
calculating an amplitude-sensitive stability boundary, such as the one 
given in Fig. 18. What is desired, then, is an approximate method that 
is capable of displaying analytically the influence of such parameters 
as ^ or Mh/a on this stability boundary. The perturbation method 
of Kamel (Ref, 1^), which was independently developed at Stanford con- 
currently with the present work, appeared to be a likely candidate. A / 
brief description of this method follows; the reader is referred to 
Ref, l4 for more details. 

Consider a given system of ordinary differential equations written 
in the following form: 


c(t;c) 


T ; € ) 


(c,t) 


k=0 


(6.1) 


Here g is a small parameter, and g is assumed to depend analytically 
on G , so that the Taylor-series expansion of Eq, (6,1) is valid. The 
idea underlying Kamel_^s method is to transform the dependent variable 
c(t) to a new set, c(t) , such that the transformed equations governing 
this new set have certain characteristics. For example, to investigate 
stability the analyst is interested in the asymptotic behavior of ?(t) , 
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He would therefore construct the transformation so that the new dependent 
variables display only the secular, or long-period, behavior of the sys- 
tem. By doing so, the' analyst presumably also comes up with a much simp- 
ler set of equations that he can solve analytically. This transformation 
is of the form 


c(T;e) = ’?(•') 


k=l 




such that c(t) satisfies 




k=0 


for the initial conditions 


( 6 . 2 ) 


(6.3) 




c(t; 0) = c' 


(6.4) 


The functions are to be determined from the requirement that the 


transformed differential equations 


c 


(.;e) . 


(6.5) 


k=0 


contain only secular or long-period terms. The ? are obtained bv 

eliminating these same terms from the governing dif fetential equations 

for the : 

k 


^ . t - S<‘) . 




i: 

i=l ' 


~ ^f,k-£ ) ^ (6*6) 
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where 


14? 









k: (i-k): 




(6.7) 


Equations (6.6) and (6.7) form recursive relations for obtaining 
and , since all other functions in these equations can be constructed 

from known quantities . This ability to obtain higher-prder approximations 
recursively is one of the principal advantages of Kamel’s method. 

With ^ and f^ calculated to any order desired, the original de- 
pendent variables c can be found in terms of c in the form 


c 






A recursive formula for the c, is 

k 


( 6 . 8 ) 


with 

\(c) = - ^(c) - 

k 

i: 

(6.9) 
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(continued) 
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and 


5c 


? - 
0,k ~ "" 


( 6 . 10 ) 


In smnrnary, then^ Kamel's algorithm proceeds by determining the 
^ ^ from Eqs , ( 6 , 6 ) and (6.7), to whatever order is desired. 

The ^ are chosen so as to cancel secular terms from Eqs. ( 6 . 6 ). These 
, in turn, give the differential equations for the transformed depen- 
dent variables, from Eqs. ( 6 . 5 ). These equations are integrated, and 
Eqs. (6.8)-(6.10) are used to determine the long-period behavior of the 
system. This method is similar in spirit to the method of multiple time 
scales, as applied by Morino (Ref. I 5 ) to panel flutter. 


6.2 Example — Two Modes 

To illustrate the application of Kamel's method, let us consider a 
two-mode version of Eqs. (2.12), with = 0 , ^ = 0 , and the two non- 
linear aerodynamic terms that were previously found to be the most impor- 
tant - those proportional to (^w^^x) and (5w/^x) (^w/^t ) . In addi- 
tion, let T = t(X )2 ^ ^ ^ ^ define 



Here allowance is being made for small perturbations from X = \ and 
from = 0 so that is the critical value of X for the linear 

theory with no in-plane loading. If substantial values of R are re- 
quired, then can be treated as another fixed system parameter, 

and terms involving it can be left with the other linear terms. With 
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these scalings, Eqs , (2.12) can be written for two inodes as 



= Ql + 


( 6 . 12 ) 
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= ^2 ^ "2 


Note that the dot superscript here denotes differentiation with respect 
to T rather than t ; this notation is restricted to this Section. 

The equations have been arranged so that all nonlinear terms are on 
the right-hand sides. With the presumption that a = 0(e., ) , these 

terms are arranged as quadratic (Q^) and cubic (C^) terms. For 


31 



[j/m = 0.01 , the characteristic equation of the linear system produces 

=1 274.5^,5 , and four roots cr^ and eigenvectors ^ 2 ^ ^ 

arranged in complex conjugate pairs as follows: 

cr^ = 1.756608 i 

= - 0.10 + 1.756606 i 

(j = - 1.756608 i 

5 

(J^ =■ - 0.10 > 1.756606 i 
U12 = - 0.9978764 + 0.0651227 i 

^22 " ■ 0.9978751 - 0.0651227 i 

" ” 0.9978764 - 0.0651227 i 
^2 " “ 0 . 9978751 + 0.0651227 i 

Since the nonlinear terms in Eqs . ( 6 , 12 ) contain time derivatives, 
the equations must be manipulated so as to put them in a form suitable 
for the application of Kamel's method. This is accomplished by vari- 

(6.15) 

( 6 . 14 ) 


ation of parameters. With ^ = (a^ ^ ^ 
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It is then required that 
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which implies 
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Substituting for the in terms of the in Eqs , (6,12) yields 
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k=l 


cr. T 
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^k = ^ = (Qi + Cl > Q2 + ^2) 


(6,16) 


Equations ( 6 , 15 ) and ( 6 ,l 6 ) constitute four equations in the four new 

* 

dependent variables 5 note that c^ does not appear in either 

the or the ^ since from Eqs. (6.l4) the a^ depend only on 

the c . Equations ( 6 .I 5 ) and ( 6 .I 6 ) can now be inverted to give a 

K 

set of equations in the desired form: 
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The matrix [B] is a 4 x 4 matrix whose inverse is 


[B] 


-1 


1 

1 

1 

1 

"12 

^22 

^52 

\2 

<^1 

^2 




°^2^22 




( 6 . 18 ) 


and [g] is a 4 x 1 column matrix 


{G} 


0 Qi + Cj^ Qg 


+ c. 


( 6 . 19 ) 
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For example, 



Here is used to denote that portion of that is quadratic in the 

a^ and therefore quadratic in the • 

. Clearly, Eqs , (6.17) cannot be integrated in closed form, although 
presumably they could be integrated numerically. A transformed set of 
dependent variables is now sought that will exhibit the asymptotic be- 
havior of the panel: 



^l,k 


( 6 . 21 ) 


The first step is to go to Eqs. (6.6) for k = 1 .A typical equation 
that results is 


1^ 



( 6 . 22 ) 


where gj ^ is now written in terms of the rather than the c 

A k k * 

The function f^^^^ is to be chosen so as to eliminate secular terms from 
Eq. (6.22). Examination of Eq. (6.20) shows that such terms will come 
about only where k = 1 in the two summations, as a result of cancel- 
lation of the exponentials. With the values of the o, being either 
pure imaginary or complex with negative real parts, there are no other 
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terms in Eq. (6.20) that could produce asymptotic instability. Hence, 
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and similar reasoning for gives 
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For k = 1 , Eqs . (6.6) then become 
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These equations are integrated to give the W for input to Eqs. (6.6) 

LfZ 

with k = 2 , and and f are chosen to eliminate any secular 

21 22 

terms. Cubic terms (s^ §2 ^ involved, and the algebra 

becomes formidable, so it will not be repeated here. (it was, in fact, 
performed on the computer.) The type of secular term that is found is 

illustrated by the form chosen for f^., ; 

.21 
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(6.26) 


The parameter is a complex constant depending in a complicated man- 

ner upon the system parameters and the elements of the matrix [B] 

For the panel problem, this second pass through the algorithm is all 
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that is needed, and the differential equation for can be written as 



^11 ^21 


Vl + 


(6.27) 


Note that is used to represent the coefficient of given in 

Eq. (6,25). Equation (6.27) has the solution 


= c 


w P " “) 

L IR -* 


. 19^ 


(6.28) 


where is the real part of and ^ is the initial magnitude 

of f |c^(o)| . The phase angle 9 ^(t) is given by = ^^j^(O)) 


^1 ” ^ 1,0 


21 


c 

2R 1,0 


In 


1 - 




•’iR ' ' 


(6.29) 


It can already be seen that one important goal has been reached, in that 
the dependence of the asymptotic behavior of the panel on the initial 
conditions is explicitly displayed in Eq. (6,28). 

Equations (6,8)-(6,10) are now utilized to construct the y 

from which the found from Eqs . (6.I5). For the c^ , we 

find 

<=kP := - ”2,k(v") 

aw k(v^) _ 

- ~ Wj k(‘=fc^'f) > (6.50) 

6 *1. 

It is not necessary, however, to calculate all of the terms in these 
equations, if stability is of interest. All that is necessary is to 
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examine the secular terms, which are contained in ♦ This fact, 

coupled with the requirement that the transformations of Kqs , ( 6 .I 5 ). 
produce real from complex ^ ^ and , permits the asymp- 

totic behavior of the system to be examined very easily. For example, 
the asymptotic behavior of a^ is given by 





(6.51) 


Using Eq. ( 6 . 28 ) gives 




-C - e 
1,0 


2R 1, 




with ~ ^1 * This equation can be manipulated to give 



where h - T) /d , -phe initial amplitude c. is seen to be one 
IK £-K 

half of ‘ The parameter is linear in and 6^ , while 

does not depend on or , so it can easily be shown that 

can be expanded in a power series in e_ , as was assumed initially 

is assigned 0(e^)) . The parameter Mh/a(o:)^ , which is the prin- 

cipal parameter governing the amplitude -sensitive instability, influences 
only Dg . ; ' 

The stability characteristics can be evaluated by examining the 
amplitude term in Eq. (6.55)« A number of cases can be distinguished: 


3T 



I . \ < 0 


For T = 0 , both numerator and denominator are negative, so 

the ratio is positive. If c. = - D , then a limit cycle 

1^0 R 

of constant amplitude 2c is obtained, regardless of the 

JL ^ L/ 

sign of D 


IR ’ 
^IR ° 


“ %o < I^rI 


, there is a limit cycle whose ampli 


tude approaches 2c. as t approaches infinity. 

“2 I I 

If c- > I D I , then there is a divergent insta- 

J. ^ vj R ' 

bility at a finite time 


In I 1 + 


h.J 


D 


2D 


IR 


B. : 


(There is a limit cycle whose amplitude approaches 
(- , independently of the magnitude of cf . 


II. > 0 : 

K 


For T = 0 , both numerator and denominator are positive. 

A. < 0 : 

Stability is indicated, since 0 

B. D^^>0 : 

There is a divergent instability at 

In / 1 + 


A_\ 


D 


2D 


IR 
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Note that Dj^ will change sign^ since it is proportional to , which 

changes sign as X goes through * 

For a numerical example, let R = 0 , so that = - e. , With 

1- X ■ . 2 1 

Mh/a(cc)^ = 0.158 , we find 


= (11.76 - 0.5500 i) , ^2 " ^ ' 1.^0^ i 


For X < X^ , 6 ^ < 0 , so this is case IIA. The panel is stable; no 
amplitude-sensitive instability can occur. The sign of indi- 
cates an increase in frequency. For X > X^ , > 0 , and case IB is 

applicable -T there is a limit cycle whose amplitude is independent pf 


c- , and the frequency is reduced from its value at X - X . (] 

that the special case c. _ = |h„j does not alter this conclusion.) 

1 ^ U K ' 

Although this numerical example did not produce an amplitude- 


sensitive instability, there are parameter regions where such an insta- 


bility will be found. The conclusions reached above are consistent with 
physical intuition, in that the changes in frequency and the character 
of the supercritical limit cycle are what would be expected in the absence 
of amplitude-sensitive instability. 
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VII . CONCLUDING REMAJUCS 


It has been shown theoretically that, under proper conditions, non- 
linear aerodynamic loading can produce unstable panel motion in a param- 
eter region that would be a stable one on the basis of a model with lin- 
ear aerodynamics . The principal factors in determining the likelihood 
of such an instability are the excitation level that the panel is ex- 
pected to encounter and the importance of the nonlinear aerodynamic 
loading in comparison with the stabilizing effect of in-plane stretching 
in the panel. The^ latter factor is measured directly by the interaction 
parameter Mh/a(o^)^ . For a given initial excitation, the critical 

value of I varies quite smoothly with this parameter. On the other 
hand, the dependence of the instability on the nature of the initial 
conditions is quite complicated. This particular conclusion should not 
be too surprising, since the system itself is far from simple. 

To the best of the author’s knowledge, amplitude -sensitive flutter 
has not been observed experimentally, at least for flat or slightly 
curved panels. Behavior of this sort is consistent with that noted in 
the experiments of Ref. 1. However, attempts to reproduce this behavior 
theoretically, with system parameters based on the experiments, were not 
successful. 

Since the type of panel considered here is far from a realistic one, 
no firm conclusions can be drawn regarding the practical impact of non- 
linear hypersonic loading. However, it is clear that any assessment of 
this effect should include an accurate determination of the true in-plane 
restraint condition. A flat panel of finite span, which is also stabilized 
by stretching in the spanwise direction, should be less susceptible to non- 
linear hypersonic effects in comparison with its two-dimensional counter- 
part. The introduction of curvature, particularly as far as shells are 
concerned, presents a different situation. It has been found (Ref. l6) 
that a nonlinear structural shell model, with linear aerodynamic loading, 
exhibits a "softening" nonlinear behavior of the same type as is produced 
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by nonlinear aerodynamic loading for the flat panels studied in this 
paper. The introduction of nonlinear aerodynamic loading into a shell 
analysis could very well reinforce significantly this behavior. Finally, 
there is the question of viscous aerodynamic effects. Within a linear 
framework, viscous effects on panel flutter are generally thought to be 
most important for Mach numbers near unity (Ref. 17). It is simply not 
known if this situation will change when nonlinear aerodynamic loading 
is significant. 

A new perturbation method for nonlinear oscillations has been ap- 
plied to a two -mode model of a panel with hypersonic aerodynamic loading. 
This method is capable of generating recursively approximations to the 
as 5 nnp tot ic behavior of the panel. The differential equations produced 
are simple enough to be integrated analytically, although the algebra 
involved is quite complicated. Fortunately, there is enough organiza- 
tion to the algebra that it can be performed relatively easily by com- 
puters, Once this software is set up, parameter surveys can readily be 
performed, and the variation of the panel stability characteristic with 
system parameters is displayed with enhanced clarity in comparison with 
the information obtained from direct integration. Although the analyti- 
cal results are capable of demonstrating the amplitude-sensitive insta- 
bility, the one numerical example that was attempted indicated that for 
the parameters chosen this instability would not occur, A final judgment 
concerning the utility and accuracy of the perturbation method must await 
detailed comparisons between the new method's predictions and those ob- 
tained from direct integration. 
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Figure 2. Dimensionless panel displacement at x = 0.75a vs. 'dimen- 
sionless time t; A = AOO, ^ = 0.01, a s: 1.0, N = 6, 
a^(0) = 0,1, linear aerodynamic terms. 
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Figure 3. Continuation; response of panel of Fig. 2. 


1.2bU 



Figure Dimensionless panel displacement at x s 0,75a vs. dimen- 
sionless time T for zero system damping; A = AOO, 

~ = 0,0^ a - 1,0, N = 6, a^(0) a 0,1, linear aerody- 
namic terms . 
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Figure J. Dimensionless panel displacement at x = 0,75a vs. dimen- 
sionless time T for zero system darrping; A t: 346. 1^ 

^ = 0,0;i a = 1,0, N = 2, a^(0) = - ^^(O) = 0.363; linear 
aerodynamic terms , 
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Figure 8, 


Dimensionless panel displacement at x = 0.75a vs. dimen- 
sionless time T for zero system damping and subcritical 
A; A = 330^ ^ = 0.0, a = 1.0, N = 6, a^(0) = - a2(0) 
=s 0.1, linear aerodynamic terms. 
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Figure 9. Response of panel of Pig. 8 -with different initial condi- 
tions — aj^(O) = - ^^{ 0 ) = 0 . 5 . 







Figure 10, Dimensionless panel displacement at x = O.J^a. vs, dimen 
sionless time t; A = 330, ^ = 0.01, Oi ~ 0.0, N - 6, 
^ = 0 . 05 , a^(0) = - ^^(O) = 1 . 71 ^ linear aerodynamic 

terms plus term proportional to (^) from (F„), , 
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Figure 12. Maximum absolute value of dimensionless panel displacement 
w at x/a = 0.75 versus dimensionless dynamic pressure 
N = 6, = 0, Oi = 1, ft/M * 0.01, Ap = 0, Mh/a - 0.05. 

Nonlinear aerodynamic loading made up of piston- theory 
terms proportional to (Sw/dx)^ and (cfw/dx) (5w/Bt) . 
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Figure 13, Maximum absolute value of dimensionless panel displacement 
w at x/a = 0,75 versus in-plane restraint parameter O. 
for X = 550. Other parameters and curve legend same as 
in Fig. 12, 
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Dimensionless stress (j at z - ±h/2 versus chord distance 
x/a, for Oi =0.1, \ ~ 550, and other parameters as listed 
for Fig. 12. Stress distributions correspond to peak -displace- 
ments plotted in Fig. 13, with that due to nonlinear aerodynami 
loading calculated for negative (into-cavity) peak. 





Figure 15, Linear stability boundary for a panel (plate-column) on 

hinged supports. Ordinate is dimensionless dynamic pressure, 
abscissa is dimensionless in-plane applied load (negative 
when panel is in compression). (After Ref. 11.) 





Figure 16, Dimensionless panel displacement w at x/a == 0.75, versus 
dimensionless time T, for point A in Fig. 15 after 
increased from point B; N ~ 6, \ = 260, a - 0,1, /x/M = 0.01, 
Ap = 0, Mh/a = 0.05, nonlinear aerodynamic loading given by 
piston-theory terms proportional to (Sfw/cx)^ and (2^/Sx) 
(Sw/St) . 







Figure 18. Comparison of linear and nonlinear stability boundaries. N = 6, ia/M = 0. 01, Mh/a = 0. 05, 
Of = 0. 1, = 0, = 1750 with a^|^(0) = -ag^O) 0, nonlinear aerodynamic terms 

proportional to (cfw/Ox)^ and (9w/0x)(8w/9t). 
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Figure 19. Variation of critical value of X with initial energy and modal content for = 0; other 
parameters same as in Fig. 18, 
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Figure 20. Variation of critical value of X with dimensionless static 
pressure difference for fixed initial conditions, R =0; 
other parameters same as in Fig. 18. 
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Figure 21. Variation of critical value of X with nonlinear interaction 
parameter for E = 1750, other parameters same as in 
Fig.. 18. ° 
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